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ABSTRACT 


Past research by Harrison and Johnston on the stability 
of pipe flow yielded only tenuous results owing to errors in 
setup of the problem and in formulation of the complex axis 
boundary conditions. 

Recent advances in the formulation of these boundary 
conditions and application of generalized stability criteria 
allowed an accurate numerical solution to be made for angular 
wave number zero. The results show that flow for this case 
is characterized by certain instabilities that have not been 
previously identified in linearized studies of this type. 

A nonuniform computational mesh was developed which 
provided dramatic reductions in computational time on a 
limited basis. 

Two data reduction programs were also developed to process 


and display data generated by the main program. 
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TABLE OF SYMBOLS 
Constant in non-uniform mesh functions given 
by equations (C-32) and (C-40) 
Partial derivatives with respect to r. 
Partial derivatives with respect to n. 
Base of natural logarithms. 


Unit vectors along the x, r and 6 axes in 
Sroadrieal coordinates. 


Cenpemenes OL the velocity vector potential 
defined in equation (2-6). 


* * 
COGELICLeEnts Of DO, D 20, is Ia, eQuatiens )(C—9 ) 
through (C-12) as defined in equations (C-13) 
ieme@ueyne (G12 2 ) 


+/-1, the imaginary unit. Also used as an index 
im Section III and Appendix D. 


iMicmiGmeer OL IMcerlor pornts in. the finite 
@iirerence mesh Of Section mil. 


Symbol denoting the phrase "of order". 
The component of the velocity vector potential 
derived from the component H by the change of 


varlable, H = xrQ. 


Reynolds number based on mean velocity and pipe 
radius. 


mae 


The streamwise velocity in Pipe Poiseuille Flow 
as defined by equation (2-11). 


Components of the complex perturbation velocity 
defined in equation (D-1). 


Complex vector potential of perturbation velocity 
defined in equation (D-2). 


G@eaindrical coordinates. 


a. tia_. Complex wave number of the perturbation 
Piero madi rect1on. 





in. Complex wave number of the perturbation 
MmehneoweEct Gece Gn where mm = 0,1,2,3, 


Ist) ene tor To increment. in the finite 
difference approximations of the derivatives 
Ot uOx 


The independent variable replacing r in the 
nonuniform mesh of Appendix C. 


Ye + 1Yy- Complex frequency of the perturbation. 
The vorticity transport equation expressed 
in abbreviated notation as defined in equation 


(ZEMEe 


The components of TF in cylindrical coordinates 
as defined in equation (2-7). 


Mesh offset parameter as defined in equations 
(G€-32) and (C-40). 


PineasmeyecCtOnL ODerarenr s(napia) 
Veeror GlLOSS-PrOduce, GCperacor. 
Brackets enclosing a matrix. 


Brackets enclosing a cOolumm vector. 





i; INTRODUCTION 


tres proolem Of finding an analytical solution to the 
pipe flow stability problem has been pursued actively ever 
Since the classical experiments of Osborne Reynolds [10] 
about 100 years ago. Up to now, however, no investigation 
Meaeebecn able to satisfactorily predict flow instabilities, 
although many approaches have been taken. 

Ssalwen and Grosch [ll] studied pipe flow with various 
angular wave numbers and sinusoidal streamwise perturba- 
tions and concluded that it was stable for all axial and 
angular wave numbers. Perturbations with exponential 
growth in space but a purely sinusoidal time variation were 
researched by Garg and Rouleau [2] and those with both 
exponential growth in space and in time by Gill [3]. Both 
concluded that the flows were stable. 

Beeawse Of EhiS inability of Linear theory to account 
for experimental fact, explanations by Davey and Drazin [1] 
involving finite disturbances and by Huang and Chen [5] 
and Leite [7] involving conditions at the pipe entrance have 
been offered. While these investigations have indeed shown 
instabilities to exist, a completely general solution to 
the linear problem has never been achieved. | 

Recently a more general theory was presented by Harrison 
[4] and further investigated by Jonnston [6]. These two 


Seueies, however, failed tO produce conclusive results due 





to mathematical errors in the problem setup and inadequate 
Pomemulatton Of the boundary conditions at the axis. Gawain 
[9] has subsequently formulated the axis boundary conditions 
in a new way which corrects the previous discrepancies and 
promises further advances. 

Por angular wave number, n, equal to zero, radical 
Simplifications result in the governing equations (Section 
Hi), indicating that this case should be approached first. 
This investigation centers on that case. 

Preliminary checks using the computer program of Ref. 6 
revealed that, of the two eigenfunctions, G and H, which 
occur in this problem and which are uncoupled for n = 0, the 
latter appeared to be the more critical. Hence the present 
research was arbitrarily restricted to investigation of the 


aeagltey Of Clgenfunction H. A similar study of the other 


eigenfunction, G, for n = 0 remains to be completed at some 
future time. Comparable calculations for other wave numbers 
Ooe— 2,3; ..-) also remain to be accomplished in the future. 


Extensive and systematic calculations of this type will be 
essential to provide the factual basis for a comprehensive 
theory of pipe flow stability. 

Reverting to the case at hand, eigenfunction H for wave 
number n = 0, we note that the program of Ref. 6 was rewritten 
er this case, incorporating the newly formulated boundary 
conditions of Ref. 9. In addition, a new, generalized 
Stability criteria was adopted. Moreover, a new technique 
was introduced which allows the use of nonuniform meshes 


to reduce computational time. 





Lastly, two data reduction programs were written to 


process data produced by the main investigative program. 


ial 





Pe lee VeonRmleriy TRANSPORT EOUATION 


Although a complete treatment of this subject is con- 
tained in Appendix A of Ref. 4 and further addressed in 
Ref. 6 and Ref. 9, it is felt that a brief overview is still 
required here to maintain continuity with previously 
referenced works. This discussion 1s an abbreviated 
memsten Of Section II of Ref. 6. 

Laminar flow of an incompressible fluid of constant 
viscosity is governed by the Navier-Stokes equation and 
the continuity equation. Taking the curl (V*) of the 
Navier-Stokes equation and introducing a perturbation 
wedioctty (v7) and vorticity (w) gives the vorticity trans- 
port equation which is equation (A-10) of Appendix A, 
Ref. 4. 

Peeeiessing this equation in terms of the complex 


velocity vector potential, W, gives 


W(x,r,8,t) = (@,F(r) + 2@,G(r) + @,H(r))en (2-1) 
where 
KX = ax + 66 + Yt iZ=72)) 
and 


v= VxwW (2-3) 





oe =) veo oe . Ley 


It should also be noted that, as shown in part one of 
Appendix G in Ref. 4, a and y are complex while 8 is a 


purely imaginary quantity defined by 


6 = en Biel 25 Sener. (2-5) 


When expressed in the form of equation (2-1), the 
vorticity transport equation becomes three simultaneous 


fourth-order differential equations of the form 


1B) ie Bs 2 


[M DG + [M. ] D~G + [M 


ID) Tal IDS tal 2 


DF |F Dir 
DG + [Mo] G = Y C{IN,] Die 
Bis H IB) Tiel 

DF F 0 

DH H 0 (2-6) 
Equations (2-5) may be further expressed in the 


abbreviated form 


eS 





Xx 

T = if = = 
ie 0 2c 1a) 
ly 0 


where [ appears to be a set of three coupled equations in 
the components of W. As eigen in ADPendd aaamom Retr a 64), 
equations (2-7) actually represent only two independent 
conditions and by an appropriate linear combination of 
is and gs equations (2-6) can be expressed as a set of 
two equations in three unknowns. The appropriate linear 


combination is given in Appendix B of Ref. 4 and yields 


the set of equations 


(2=18)) 


Except for the case where n is equal to zero, equations 
(2-8) do not uncouple. The linear combination given by the 
second of equations (2-8) does, however, reduce the highest 
Speer derivative of G{r) in equations (2-6) to second order. 
Appendix C of Ref. 4 illustrates the redundancy of the three 
components of W, allowing one of these components to be 
arbitrarily set to zero for all r. The maximum benefits 


of equations (2-8) are obtained if 


etic) = 0 ea) 





Incorporating equations (2-8) and (2-9) into equations 


We) results in the form 


“ 


Dire p°c D°G 
[Mt] ee Pees | fo (MoM! | 
2 DH 3 D3H 7 DH 
DG G D*c 
ao (M | + PME = Sp Oa 
1 DH 2 H 2 D-H 
DG G 0 
et [Ni ] ate [N‘ ] ) = (2=10) 
l DH H 0 


where the coefficient matrices are given by equations (2-10) 
meeemema ({2—-1/) of Ref. 6. It 18S appropriate to note that 
these same coefficient matrices appear in Ref. 9, equations 
Peesthroucgh (A9), in a slightly different form resulting 


from the substitutions 


Qe Seba eee Caine 
2. B° 
VE = eo! ate 7) and (2-12) 
te 
io) ee 
T = oU- =(a* + By) . (aes 
e 1g 


As discussed in the previous section, the case where 


82 = in, ioe — 18, (2~-14) 


15 





Heads tO great Simplifications in equations (2-10), 

(2-12) and (2-13). In particular, equations (2-10) 

uncouple and allow an independent investigation of either 
[meee ewor = AS a result of the findings discussed in Section I, 
it was decided to explore the function H only. This reduced 
equation (2-10) to that of equation (A-6) of Appendix A, 
which is a linear, homogeneous fourth order differential 


Sailatien in H(r). 


16 





III. NUMERICAL METHODS 


Substituting the change of variable H = rQ as given 
in equation (A-1) and the coefficients defined in equations 
(A-11) through (A-18) into the vorticity transport relation, 


equation (A-6), gives the expression 
M ae Ml p39 ae eal D“0 em eo) yeeros 
4 3 2 1 0 
- y[N,D°Q + N,DQ +N,Q] = 0 321 
2 1 0 ' oe 
which is a homogenecus fourth order differential equation 


mamere). The boundary conditions for this case are 


fetmivead in detail in Ref. 9 as 


ii == 0 
DOK) ee 
(3-2) 
DQ(0) = 0 
3) 0) (0) ae 


The boundary finite difference equations derived in Appendix 
Beecremecauations (3-2), along with the standard central 
difference equations given in Ref. 6, allow the function Q(r) 
to be approximated by a finite number of discrete unknowns. 
As shown by Figure 3-1 below, the non-dimensionalized 


radius of the pipe is divided into a one-dimensional 


sy 





r=0Oe¢eeeveeereeveeeeeeeeeneeeee ee eee et 


@ | 
@ 2 


@ N- 
@ N 
=D COCCMWPDWG GG  eeeWe.§. 


Figure 3-1 Finite Difference Mesh 


computational mesh consisting of N interior points, N+tl 
Mieeeyals, and N+2 total points, including the boundary 
Points at tr = 1 and r= 0. As will be discussed later, 
the spacing between these points may or may not be uniform. 


For the uniform case, the spacing is defined by 


6 l=, i) ee (3=3) 


For the nonuniform case, a change of independent variable 
is performed. The spacing of the new independent variable, 
fpeesesti li given by equation (3-3). 

With a nonuniform mesh, the points shown in Figure 3-1 


Will be concentrated near the axis or near the wall according 





tO the type of offset specified. These effects are dis- 
cussed in detail in Section IV. 

Substitution of the finite difference equations of 
Appendix B into equation (3-1) results in a set of N, 
linear, algebraic difference equations in terms of the 
unknown value of Q at each of the N interior points of 
the computational mesh. Since each of these equations 
Moore cine Gem Of a linear combination of the ith, central, 
point and the two, three or four adjacent points (depending 
on the order of the derivative being approximated), this 
system of equations consists of a coefficient array multi- 
plying a vector containing the unknown value of the function 
Q at each of the N interior points. This technique allows 
the problem to be converted into an eigenvalue problem of 


Eae f£orm 


Rote ¥ [10 ae a0 (3-4) 


Paeweeae basic composition of the arrays [X] and [Y] and 
the vector {Q} as illustrated in Figure 3-2 below. 

It should be noted at this point that Figure 3-2 differs 
Somewhat from the normal finite difference banded matrix in 
the first two rows and last row because of the method of 
deriving the finite difference approximations at the boun- 
daries. Additionally, the order of the N unknowns has been 
reversed from that of Ref. 6. This was done to conform to 


Seandard matrix notation. 








e©e2ee8ee 

eeeoe5e Q, 
eee ee QO; 
eececee 
ae on ‘ at 
eee e @ ; 
eee e : 

e @ @ QO, -; 
ee 0 Qn 





maguire 3—2jbas lewecompos ition. ofeCoefficient 
AEeays Ange ecteougeol Unknowns 


This array is established by the subroutine MSET2 in 
conjunction with the subroutine MSET1 and function sub- 
programs CQMIE1 and CQM2E1, which compute the numerical 
value for each element in the array. Subroutine MSETI pro- 
vides the coefficients given by equations (A-11) through 
eeoveor Appendix A or by equations (C-24) through (C-31) 
mena mOnuUni form mesh is specified. Function CQMIE] then 
computes the values for each of the elements of array [X] 
in equation (3-4) using the coefficients passed from sub- 
meomEadie Monrl in vector COM1. Function subprogram COQM2E1 
Peeeemms the same function for matrix [Y] in equation (3-4) 
uSing the coefficients passed in vector CQM2. 

The solution of the eigenvalue problem as formulated 
~emehesm@Gine 15 Carried out by the controlling subroutine 
of program PIPEOQ, subroutine STAB, by the following steps: 

1) Subroutine MSET2 is called twice to set up the coeffi- 


cient matrices [X] and [Y] of equation (3-4). 


AK 





Subroutine CDMTIN is then called to invert matrix 
[Y], the second coefficient array in equation (3-4). 
CDMTIN was obtained from the IBM Library routine 
CMTRIN by modifying it to accept double precision 
arrays. 

Both coefficient arrays, [X] and [Y], are then pee= 
motlet>iied by [Y] -. Since multiplication of an 
array by its inverse invariably results in the 


iemeiay Matrix, (1), only the DEeauct iy] - 


xs 
computed using subroutine MULM. This converts the 
eigenvalue problem of equation (3-4) to the more 


Sea eneronal form 

eZ) “SaelliyiO) == 0 (Sa 5) 
where 
[X] C36) 


Since all programs currently available for solving 
equations (3-5) require that the real and imaginary 
parts of the elements of [Z] be presented in separate 
arrays, subroutine DSPLIT is called to accomplish 
elolakss © 

The eigenvalues and eigenvectors of equations (3-5) 
SGeneempucted Using subroutines EBALAC, EHESSC, ELRH2C 


and EBBCKC which are available through the International 
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Meieeane ostabistics tibrary. Subroutine BBALAC 
Paeamees Matrix [4] by@equalizing the exponents of 
all terms. The details of this transformation are 
retained for later use. The balanced matrix is then 
passed to subroutine EHESSC where it is reduced 
into the complex upper Hessenberg form. Subroutine 
ELRH2C then solves for the eigenvalues and eigen- 
Vectors. To’ Eransiorm the elgenvectors back into 
the original unbalanced form, EBBCKC is finally 
called using information passed from subroutine 
BBALAC. 

For each solution, subroutine STAB determines the 
least stable eigenvalue (largest algebraic value) and 
Ree eee phe aT 
and KSET to file FTO2F001. The eigenvector corres- 


then writes the values of N, Ror ei 


ponding to the least stable eigenvalue is also written 
to FILE FTO02F001 when MODENO is set equal to one. 
Control of subroutine STAB is accomplished by the 
main program, PIPEO. This program is a time-sharing 
(CP/CMS) program. Modes one and three compute the 
Stability of the flow for a given set of input con- 
ditions. Mode one writes the least stable eigenvector 
to FILE FTO2FO001 while this output is inhibited when 
MODENO is set equal to three. To generate data for 
pecqmeanenliGheN, program PIPEO Must Be run with 


MODENO equal to one. 


Ze 





Mode two operation generates a grid of stability 
values (stability map) based on parameters read in 
from FILE FTOIFO0O1. Due to the long run time in 
this mode, only small meshes can be generated under 
CP/CMS. Longer ee must be accomplished under batch, 
with changes to the program as specified in the comments 
section. Data is output to file FTO3F001 when MODENO 
1s equal to two and ts comvatible with program STBCONT. 
The plotting programs EIGFCN and STBCONT were used 
to process the data generated by program PIPEO in 
modes one and three, respectively. Program EIGFCN 
generates normalized plots of the perturbation velocity, 
Maraoewa Lunction Of radius, rr. The perturbation 
velocities generated in accordance with Appendix D 
were normalized in two steps. First the perturbation 
velocity of largest magnitude was determined. Letting 
this velocity be termed Yas a normalizing constant 
producing unit magnitude and zero phase angle in U, 


Wesefouna i2n the following manner: 


ies 
Ce om u + iu. ; (3-7) 


then 


Cu. = It se 10) (Ses) 


ae 





a) 


where C is the normalizing constant. Thus, 





me aE - sislel @ tye 
u te Jae 7 ? (see!) 
RC ui Ge (u - u.-) 
Re a1 @ 
= 
- C 
3 ie 2 (3-10) 
GC 


where By is the complex conjugate of up. 

The nondimensionalized radius values were taken 
@irectiy from the data cards for uniform meshes or 
computed from equations (C-32) or (C-40) in the case 
of a nonuniform mesh. 

Eaejmadil olBCONT plots the stability contours against 
Op and Ay: The stability map generated by program 
PIPEO is searched columnwise and rowwise for sign 
changes for each of the three stability criteria 
discussed in Section V and Ref. 9. The points are 
etemolLocted, producing contours Of incipient, criti- 
cal and fully developed instability and areas that 
denote stable flow and subcritical, supercritical and 
meaweermeritical instability. 

Both programs, EIGFCN and STBCONT, utilize the NPS 
WemRosatnhG plotter, certain Duilt-imysVERSATEC sub- 


routines, and subroutine PLOTG. These routines are 


only accessable when running under FORTCLGW. 
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IV. RESULTS 


A. STABILITY 


Since an understanding of the term stability is necessary 
fomemeerpret the results of this investigation, a brief 
discussion is presented here. A complete discussion of 
eiemegemeralized criteria of stability is given by Gawain 
[SNe 

Wiemenaracteristics of the flow for the case n = 0 are 
set by the parameters R. and a. For fixed values of these 
parameters, the solution of equations (3-5) is a set of N 
eigenvalues, y, and their corresponding eigenvectors, Q. 

As can readily be seen from equation (2-1), the value of 

the real part of the complex eigenvalue y will determine 

the growth or decay rate in time of the perturbation. Since 
positive values of the real part of y represent an exponen- 
eiategrowth rate in time, the most important y is the one 
having the largest algebraic value for its real part. This 
root is termed the least stable root and will be represented 
by the symbol Yer: As the stability represented by Yer 10) 
that seen by a fixed obServer, it is not the most general 
Gamreriton. As derived in Ref. 9, a more appropriate stability 
criterion is that based on an axis system moving at the 
average volumetric velocity of the flow. This criteria is 


termed and is defined by Ref. 9 as 


* 
RL 
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ey P (4-1) 


For this and subsequent discussions, the subscript will be 
dropped and y will refer to the quantity defined by 
equation (4-1). Three stability cases arise from this 
equation. The first is termed incipient instability and 


is defined by 
The second case, termed critical instability, is given by 
1 = (4-3) 


and, lastly, the case termed fully developed instability 


1s Said to exist when 
(4-4) 


Pie ecransition from stable Elow to fully developed 
instability is progressive and several distinct stages are 
given in Ref. 9 to describe this transition. The region 
Past@eimemolent to critical instability 1s termed subcriti- 
aeueeerstao.lity, that from critical instability to fully 
developed instability is called supercritical instability 
while that beyond fully developed instability is termed 


moereritical instability. 





Bee PERTURBATION VELOCITY PLOTS 

Initial investigation of the function Q was centered 
around plotting its appearance in the region of interest. 

A Reynolds number of 1150 (2300 based on diameter) was 
chosen as this value is generally accepted as the nominal 
Peas FOr transition to turbulent flow. The value of a was 
set eee o-> t+ 2 10.0 for the major part of the investiga-— 
€10n aS preliminary checks revealed that supercritical 
instabilities were present for this value. A secondary 
Reynolds number of 4000 was chosen to show trends. 

The quantity chosen as the most realistic and represen- 
tative of the eigenfunction Q is the axial perturbation 
velocity, u. This quantity was derived from the elements 
of the least stable eigenvector as outlined in Appendix D. 
tabelally, Re and oe were held fixed and dp was varied 
Over a range of positive and negative values. For values 
Gut Op below about two, the normalized perturbation velocity 
was found to have all activity near the axis with a decay 
essentially to zero by r = 0.3. A typical plot of u versus 


r for an a, in this range is shown in Figure 4-l. When 


R 
Ap was made sufficiently positive, the plot changed signi- 
ficantly in both appearance and region of activity. Figure 
4-2 shows a plot of u for op = Zo ene waAcelVity Can now 
be seen to be concentrated near the wall, with most of the 
qemieity occurring at ©r values greater than 0.7. 


Pienetigniene: Particular Selationshtpe becween the mature 


of u and the stability of the flow was evident or expected, 





the plots were nevertheless valuable as indicators for 
various parameters involved in the investigation. 

First, as can be seen by the differences in PaGuLes 
4-] and 4-2, the plots were ideal indicators of changes in 
the nature of the function Q. Secondly, the adequacy of 
the mesh could be directly observed by noting the number 
of points defining the curves in regions of high activity. 
Figures 4-3, 4-4 and 4-5 show the same conditions as 
Figure 4-1 but with decreasing number of mesh points, N. 
Lastly, the effects of nonuniform meshes could be observed 


aS will be discussed later in this section. 


er  SPeaGbiLiTy CONTOUR PLOTS 

The principal results of this investigation are shown 
in Figures 4-6 and 4-7. Although these two figures pertain 
to only a limited portion of the complex a plane, they do 
represent a Significant advance in the investigation of 
pipe flow stability. As can be seen in these figures, the 
Miers characterized by regions of differing stability, 
memerng £rom stable through supercritical instability. Note 
that these two figures correspond to Reynolds numbers of 
1150 and 4000, respectively. This is a result that has 
not, to this writer's knowledge, been heretofore achieved 
by a linearized analysis of fully developed pipe flow. The 
Figures also show that, as has been born out by previous 
investigations, flow for purely sinusoidal oscillations 


(O, = 0) is stable. Additionally, a comparison of Figures 4-6 
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and 4-7 shows the effect of Reynolds number on the flow 
Stability. It is clear from this comparison that an increase 
in Reynolds number reduces the size of the stable regions 

in the complex a plane; in other words, stability decreases 
with increasing Reynolds number. This trend agrees with 
meegecnieral experience pertaining to fluid flow. Lastly, 

Mee ertect of the real and imaginary parts of the wave 

number a can readily be seen. For Ops increaSingly nega- 

tive values produce successively greater levels of instability. 
While a contour plot was not produced for positive values 


St oO point checks of stability in this region suggest 


R! 
that somewhat similar contours exist in the right half- 
plane also. For are increasing values produce increasing 
Stability. This effect is also more pronounced at the lower 


Reynolds number. 


Dee NONUNTP ORM MESH EFFECTS 

One of the difficulties in this investigation was the 
relatively long computing time required to obtain an accurate 
solution, especially when operating under CP/CMS (time-sharing). 
Baewmayor factor controlling computing time was the number 
of interior mesh points, N. AS an example, an increase in 
N of 50 percent resulted in a fourfold increase in computing 
time. Therefore, the desired objectives of rapidity and 
Beeurdesy were in direct conflict. Additionally, follow- 
On investigations for values of angular wave number n other 
than zero involve matrices twice the order required for this 


case because of the coupling of equations (2-8). 





For these reasons, a nonuniform mesh was developed 
menOQeain increased accuracy at lower values of N. The 
nature of the velocities as seen in Figures 4-1 and 4-2 
shows that a high degree of resolution in the computational 
Mesh 1s only required in the vicinity of the axis (a, less 


R 


than about 2) or the wall (a, greater than about 2). It 


R 
was therefore theoretically possible to redistribute the 
points at moderate values of N to attain resolutions equiva- 
lent to much finer (and more time-consuming) uniform meshes. 
As can be seen from Figures 4-8 and 4-9, the value of 
<a varies with the number of mesh points, N. Theoretically, 
each of these curves would approach some limiting value 
if N were increased without bound, and it is this theoreti- 
ea) limit that represents the required solution. In prac- 
tice, it is adequate to approximate the unknown limit by 
a point that lies on the relatively flat portion of the 
curve at a value of N which is practically attainable and 
which does not involve a prohibitively long computing time. 
fe has been found in this investigation that N = 79 ful-, 
pies schese Conditions. 
The conversion to a nonuniform mesh involved a change 
of independent variable and the introduction of an analyti-~ 
cal function to control the distribution of the mesh points. 
The details of these steps are given in Appendix C. By 
varying the mesh offset parameter, A, it was possible to 


* ® ® 
vary ( over a wide range. To determine when the high 
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accuracy solution (N = 79) and the nonuniform solutions 
were approximately equal, Y Was plotted versus A for fixed 
values of Ra: a and N with the value of a Oi ha" 3/ 9 ais 
feeeeeorence. Figure 4-10 shows a plot of this type for 

N= 31. The appropriate value of X can be seen to be 
approximately 1.1. Figure 4-11 is the perturbation velocity 
eeeewer the Solution for N = 31 and } = 1.1 for the Same 

Re and a aS Figure 4-1. Note that the va values are equal 
for these two figures. While the resolution of Figure 4-1l 
Baenoc Guate as fine as that of Figure 4-1, a comparison 

of Figure 4-11 with Figure 4-5 makes the improved resolu- 
tion obvious. Figures 4-2 and 4-12 are similar to Figures 
4-1 and 4-ll except that a wall offset was used. Note that 
£0r this case ’} = 1.2, which points to a drawback of the 
nonuniform mesh, that of dependence on input conditions. 
While a check of A dependence on a was not made, it most 
Pm@eemly Exists. There is also, however, the possibility 
that for small regions of the complex a plane, the varia- 
tions in X are small enough to allow an average value of 

X to be nearly optimum for the entire region. While not 
used for the main results of this study, the method as 
developed here may well prove to be of maximum utility in 


follow-on investigations of higher angular wave numbers. 


Pe eNeMERICAL ACCURACY 
To ensure that the solutions presented here were of 


SuEficient accuracy, two separate checks were made. The 





ircst, Y dependence on N, is the most commonly used 
@-2Cerion. 

For a solution to be accurate, it should be Vai Cail 
independent of mesh fineness, that is, of N. The required 
magnitude of N for an accurate solution was found by plotting 
Y against N. Figures 4-8 and 4-9 both show that the 
solution is well converged for N = 79 at Reynolds numbers 
Seelt50 and 4000, as y changes by only .001 to .003 
mrom N= 31 to N = 79 for both values of Reynolds number. 

Miewsecond verification of the solution, so obvious 
that 1t 1s sometimes overlooked, involves simply substi- 
tuting the numerical solution (least stable eigenvector) 
into the governing equation to ensure that it is indeed 
being satisfied. A short program was independently written 
to check the finite difference representations of equation 
(3-1) at the first and last interior stations and at a mid- 
mBaertus Station. Initial checks of numerical solutions 
yielded unsatisfactory results and led to the discovery of 
Various programming errors. In particular, it was dis- 
covered that four double precision constants in the finite 
difference approximations were lacking the required "DO" 
exponent. Elimination of these seemingly trivial errors 
resulted in a surprising four order-of-magnitude improvement 
in the accuracy of the solution, with the left side of 


4 EOuerECe t Meo. 


equation (3-1) improving from order 10 
Mem omNSErUCtTIVe tO note at Ehis point that the order 


SEemegqnitude of the left side of equation (3-1) 1s not the 


oe 





ae 
en 


true measure of its satisfaction. A more correct procedure 
is to compare this value with the largest term in the equa- 
tion. When examined from this viewpoint, the relative 

error for solutions at ss = 1150 and iS 4000 are found to 


11 1) 10°14, a very satisfactory result. 


be of order 10° 
Therefore, by these results, the solutions presented 
here are both virtually independent of N and satisfy the 
governing differential equation to a high degree. The 
efforts expended to reach these conclusions were well 


worth the result and also point out that attention to detail 


memitanacamental to accurate numerical results. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


The implementation of the newly developed boundary 
conditions of Gawain [9] has permitted a Stable, numerical 
solution to the linearized vorticity transport equation. 
Bre results of the numerical solution are presented in 
S--210m 1V and show that the stability of pipe Poiseuille 
flow 1s governed by the three parameters, 


a ane hee 
e 


ota 
Bomeeacrelcular, both positive and negative values of ps 
that is, streamwise growth and decay in space, if suffi- 
Cclently large, produce unstable growth rates in time. This 
result is new and it is consistent with Eee Kewl se <Pcials 
Mental fact that transition to turbulent flow depends not 
only on Reynolds number but also on the general character 
Senthe perturbations which exist in the flow. 

The perturbation velocity plots of Section IV represent 
the first practical look at the function Q. These plots 
were valuable indicators for adequacy of mesh fineness, 
Piaeets, N, Changes in the nature of the function Q and 
effects of a nonuniform mesh. 

No instabilities were discovered for purely sinusoidal 
perturbations (O15 =U) hts 2S conststeeme wlen) the 
Pisevylous Investigation of Ref. ll, but should not be 


assumed for investigations of other angular wave numbers, 


ale 3. ) « 
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Adequate numerical accuracy was proven by demonstrating 
that the solution was virtually independent of the number 
of mesh points, N, and that it satisfied to a high degree 
an independent check of the governing differential equation. 
This procedure should also be carried out in future inves- 
agiatelons Prior to conducting full scale data runs. 

This study suggests that similar, and perhaps even more 
rewarding results will be obtained for the higher angular 
wave numbers. Although lengthy, programming is straightfor- 
ward 1f approached systematically. The general organization 
Of the programs of Ref. 4 or Ref. 6 should be helpful in 
this task. It is recommended that the case for n = 1 be 
undertaken as a follow-on to this study. 

The nonuniform computational mesh was shown to be a 
mewertul tool in the reduction of computational time. At 
the same time, however, the dependence of the mesh offset 
parameter, 4», On input conditions needs to be investigated 


further to realize the full potential of this technique. 





APPENDIX A 


Pent vAtLON OF VORTICITY TRANSPORT EQUATION COEFFICIENTS 


From the change of variable introduced in Ref. 9, the 


function H for the case n = 0 is expressed by 


As ee) (A-1) 
Taking derivatives 
paw =— 2DO +°O (A-2) 
D°H = rD“O + 2D0 (A-3) 
D°H = rD°0 + 3D°0 (A-4) 
Dy = rD*9 + 4D°0 (A-5) 


Let the '*' superscript denote element (2,2) of 
peices (Al) through (A9) of Ret. 9. SincesroOnr nm = 07 
only the function H was investigated, Squcacrons. (2-10) 
become 
2 


x 4 * 3 # 
[eet Me bDH+ M. DH+ M 


2S sags oe 
4 3 5 1 DH 0 =! 


~ * 
- y(N,"D°H + N, DH + No HI] = 0 ee) 


it 


Substituting for H, equation (A-6) becomes 





M "(rp494+4D39} +M 


4 * (rp70+3D°Q} + M *(rp“9+2D0} 


3 2 


¥* a 
+ M, {rDQ+Q} + Mo {rQ} - y(N,” {xD*Q+2D9} 


3 N, {xDQ+Q} - Ny (rh — an) (A-7) 


Before proceeding further, it should be noted that the 
Ref. 9 matrices from which the coefficients for equation 
(A-7) were taken were obtained from matrices (2-10) through (2-17) 


of Ref. 6 by means of the following substitutions: 


GS AO ee (A-8) 
ig 
t = ne ~~ (A-9) 
18 
gh 
1,2 2 
T = aU - gla - 5) (A-10) 
e 1g 


Defining the new coefficients for equation (A-7) as M 


0 
Earough M, and No through N. 
— = Gent 
M, = cM, = R 
e 
: i 6 (A-12) 
M. = 4M, + rM. = = a 
e 
* * ees 2: 
= = _— |} = 
M., 3M, + rM, raU Rr 2a“r} (A-13) 
* 
ye ee a es ees 
1 2 i R 2 
er 
. Oe a ar (A-15) 
My = MS + rMy =) re R 


(= 


*. 


7 


7a. 





N = rN = -—r (A-16) 


2 Z 
= D * * 
Ny = N, + rN, == 3 (A-17) 
ee = NS fae See 
9 = Ny rNy = -a r (a eaters) 


Upon making use of the foregoing substitutions, the 
governing relation can finally be reduced to the form 


Beevyvrousily shown in equation (3-1). 





REP END B 


FINITE DIFFERENCE EQUATIONS 


Improved finite difference equations for the boundaries 
were obtained by not using the virtual point method of 
Ref. 4 and Ref. 6 and deriving the forms directly from the 
boundary conditions of Appendix A. The equations thus 
formed are also of consistent order truncation error, 
Pegmbetcantly improving the accuracy of the solution [Ref. 8]. 

Because of a peculiarity in the form of the consistent 
Becena orcer truncation error equations at the axis, a 
Simgularity resulted for a equal to zero. Consistent third 
order truncation error equations eliminated this problem. 


From Appendix A, the axis boundary conditions are 


DO(0) = O aad p?0(0) = O (ea) 


Representing Q by a power series and applying equations (B-1) 


wields 
2 4 5 
- 3g 5 1% 
@( x } =e (0) + D°Q(0) 57 aD. OVO) ql Q(0)E> 
6 r° 
eee, Q(0) EF Te) ‘yore (B-2) 
Using five mesh points at r = 5, 26, 36, 46 and 56 results 


ieee matrix 











“hl a ns) 
Q,} = zag << 2 Sea) 2 O54 
4 8 SP ao | fe Pao) 
(B-3) 
Differentiating equation (B-2) and substituting r = 6 gives 
gia, Matrix form) 
Q(6) - FR Q(0) 
$DQ (5) l = Ff aoa 
oO (6) = 1 — + SB HCO) + By! 
§3p79(6) 0 1 §°D°Q (0) 
ls4p40 (8) 0 1 7 68% 0} 








Use) 


mee iAl and [B] denote the coefficient matrices of equations 


Os 
(B-3) and (B-4) respectively. The values of Q(0), S°D°Q(0), 


34p49(0), 3°p°9(0) and 6°D°Q(0) may be solved for by 





Q(0) Q 








i 
Zee 
6 “DQ (0) Q 
2 
4 4 - 
§*D*9 (0) = fa] t Q, Soy! (B-5) 
5 
§°D°Q (0) 9 
4 
66 
|: D 0(0} Q. 
Putting equation (B-5) into equation (B-4), 
6 
Q(0o) Q, 
OBO Mey), Q. 
ee 2 - 
s“p°Q(s)} = [BIA] ~ ja,f + 067 (B-6) 
$°p°9 (6) Q, 
~4 A 
$“D%Q (8) Q. 
The last line of this set of equations gives 
D0 (6) = oq (--9115646260, + 2.7502429550, = 30437317790. 
- 1.424684160, = -2196307090,) + O67 (B=7) 


To solve for DO, the rightmost column and bottom 
Soweare eCliminated from matrices [A] and [B] then these 


new matrices are inserted into equations (B-5) and (B-6). 





The bottom line of equatina (B-6) will now give the 


5 ; 
expression for D°Q(6) with a consistent third order 


EpEemmcatlon error. D“0 (6) euGie WO ey Were so lyee mom ina 


Similar manner. 


D-9(6) = 

Do = + 
6 

DQ = =(- 


1 

“3 (1. 8251655639, ~ 3.2503311260, + 1.6609271520, 
8) 

- .2357615899,) + 05? (B-8) 

385 4, 8 1 3 
( 5 Qs fe is Q. a 3p Q3) + O06 (B-9) 
2 2 3 

= 0, +£,) + 08 a0) 


Biewco the complexity of the boundary conditions, it 


mesomeeeiaedq that consistent third order truncation error 


equations should also be used at r = 26. For this the [B] 


matrix only need be changed as equation (B-2) is unchanged 


Perdis Station. The new matrix [B] is formed by differ- 


entiating equation (B-2) and making the substitution r = 26. 


Proceeding as for r = 6 gives the following finite differ- 


Suce approximations 


p*9 (26) 


p°Q(28) 


4, (-3.103401360, + 6.903012634Q,- 53422740530, 
5 

+ 1.66083577Q, ~ 0.1234207970,) + Oo meses 
=; (.8688741720, - . 9377483450. - 2543046360, 


On 


3 


cf -3231788080,) + O6 (3) 





2 eek 28 19 
Eee) “= <5 732 - 720, + 530,) + O67 (eae 
DQ(28) = F(- 40, + 40,) + 063 (B-14) 


It should aiso be noted that the value of Q at r = 0 may 


be solved for from the top line of equations (B-5) 


Q(0) = (1.795918367Q, - 1.24781341Q., + .606413994Q, 


- .177842566Q, + .023323615Q.) + os? (B-15) 


The central difference equations given by Ref. 6 were 
already consistent second order truncation error equations 
as confirmed by Ref. 8 and were retained. 

For the wall, the clamped end, consistent second order 
equations (5) through (8) of Table II, Ref. 8 were modified 
for the "right boundary" using the procedure given in Section 


5 of that reference. 


. ee g ; , 
D°Q(1-6) = =Fl- FOx_3 + FOy_2 7 9Qy-y 7 1OQy) + O6 
(B-16) 
pD°Q(1-6) = st FQu> + 304) + 06% (B-17) 
D’9(1-8) = 41, - 20.) + 08° (B~18) 
0 
Da(i-6) = S(- Z0y_,) + 08° (B-19) 





Since the wall finite difference approximations were 
@e Only Second order truncation error, the approximations 


Bor DO through re at r = 1-26 were obtained directly from 


the central difference equations with Q(1) = 0. 
PEO T= 26) = ate - 40 + 60 - 40.) + 06° (B-20) 
4 N-3 N-2 N-1 N 
3 <«, Mee! : 2 : 
Beet —26) >. = en Fy 3 +t On_5 Q,) + 06 (B-21) 
poe) = (0... - 20... +0.) + 08- (B-22) 
32 N=-2 N=-1 N 
, eee oil 2 _ 
OK 1-26) = a Fea ae aes cS (a2) 





AEE NIDUEDC 10 


NONUNIFORM MESH 


MemecOnterol Ehe distribution of a fixed number of mesh 
points, a change of the independent variable from r to n 


was performed. 


Q(n) | (one) 


1—O 
Il 


Te) Pee iy) (C-2) 


The derivative with respect to r becomes 


* is *® 
D = (Dr) *D (C-3) 
where 
bs d d 
=_- aes = — -4 
D a and D ar ‘G—4)) 
BO, p*9 ... Can now be expressed in terms of the new 
independent variable, Nn. 
* _ * 
DO = (D rr) Lp Q (C=) 
Y * -l * 
D°Q = D(DQ) = (Dr) “D (DQ) 


mee 2D “or (> A) @ “x Die (c-6) 


= 





D(D“Q) 


x 
by 1g) 


* 
ee) 


D(D70 


* 
(Dr) 


= 3 


) 


—4 


ee 
= 4D" R) 


p39 = say =e 


*4 


—4 


a 


* 
(D> 


“ 


* 
D (D“Q) 


eZ Ee 


(oes) OF ae 


2 


< =5 * 
Y= 450) (Dt) 


2 ee: 


* -5 * 
De Oe morb. 5) (Dae) Da® 


='6 


1 
eo (Dr) ( 


ei 


* 
meee tD cr) (( 


* 
yD) 6) 


ne 


-1D9 (C-7) 


D r)-4(D r) -(D -xr}]D oO 


* 
D Boa 2 


zo 


(D4r)1D0 


3 


416 (pe) i se) oy oe 


(Ga3) 


The derivatives of Q with respect to r can now be written 


where 


WEae 


Ze 


SS 


44 


* 
(D xr) 


* 
ar [oucact 


#4 
ic Oo + fy 


nae 


3) 


3 


* * 
fee) «2 Oi’ £51) fe) 


p*20 ek 


+ 
D 36 18 


ie 


42 


* 
Q 


x2 P 
iB EO as Ai Q 


(ee2) 


(G>1 87) 


(Cala) 


(Gal) 


(aus) 





f,5 = (D r)~° 

f£,, = -(Dr) 7 (D"r) 

t 33 = Op ie) 

f,, = -3(D 4)" (D'*r) 

34 = 3(D r) > (D 2r)2 = bo ay te 
in, = (OLAS 

ae = 6p =) 2s Sey 
psi) zr) (D -r)° - 4(p rc) > (D x) 
f,, = -15(D r) /(p “*r)* + 10(D r) °(D ” 


= pire) pee 


Seesetemeing equations (C-9) through (C-12) into the 


vorticity transport equation (A-6) yields 


k *4 k_ *2 . 
My IDL Se M. Da + D 


k *O 
-y [N. Be ROME a Ny 


x * 
DQ + M 


1 


k * k 
Bie (Oe is Mo @) 


xk * * 
DQrt No QO) = QO 


pay) ec (8, 


(C=14)) 


(CeiS} 


(C—16) 


re ale} 


(C=i3) 


(Gai 9)) 


(e290) 


USS Zl) 


(E=22)) 


hea 2) 





where 


M, = M,f,, (e=oa) 
* 

M, = Mt 43 = M,f., (G=25)) 
* £ 

M. = M, A? ~ M,f,, + Mof,5 (C-26) 
* 

ee oe 423, Mofo, (C= 27 | 
* 
* 

ees 22 (C-29) 
* 

Nie = Noto, * Ny f,, (C-30) 
* 

No = N Ce 2UL) 


In order to concentrate the mesh points at the axis, the 


mance On 


ie ESS ierch got nib np) ) (e=32)) 


was chosen where A iS a parameter controlling the degree 


of concentration of mesh points near the axis. Equation 


fe-—32)5must satisfy the two conditions 
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me 0 at iW =o2te (C-33) 


and 


Substituting equation (C-33) into (C-32) gives 


Gene /tanhe sk x (C-35) 


Computing derivatives 


De = Ci/coshA(1-n) : (C-36) 
Dee IC) (tanhi (1-1) /cosh-i(1-n) | (e=2m 
Me 20) -{ (1-25 inh?) (1-n) 1/cosh- i (1-n) } (C-38) 
Dr = 8Cd*{tanh°i(1-n) /cosh2) (1-n) 1 (C225) 


To shift the mesh point concentration to the wall, the 


muMectLOn 


r = Ctanh An {C-40) 


Mas selected. Satisfying equations (C=-33) and (C-34) for 


this equation also gives equation (C-35). The derivatives 





memic-40) are given by equations (C-36) through (C-39) if 
n is substituted for all occurrences of (l-n) and the signs 
of equations (C-37) and (C-39) are reversed. Figures C-l 
and C-2 show equations (C-32) and (C-40) for four selected 


values of the parameter i. 








FIGURE C-l. R Versus n for Four Selected Values 
of Lambda - Axis Offset 


6s 





Os 
Q.6 6.0 
Q.4 
Ole 2 
G0 ; 
0.0 O22 0.4 0.6 0.8 
mae 


FIGURE C-2. R Versus n for Four Selected 
Values of Lambda - Wall Offset 
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APPENDIX D 


DERIVATION OF PERTURBATION VELOCITIES 


From Ref. 4, Appendix E, equations E-6 through E-8: 


Coa) 
v(x) = {(AJw + [B]DW (D-1) 
w (5) 
eae i 
0 = = F 0 0 ab [ee 
= = ie G ot) 0 0 DG 
0 x 0 | | “Si prey DG| 
Ce), 
For this case 8 = ni = 0 and F = DF = 0. Restricting the 


investigation to the function H for the reason expressed 


mimcection £ and solving for u(r) gives 


(ig) Se EL (D2) 


Performing the change of variable 


H = rQ (Das) 


DH = 0 + FDO (D-5) 





Ee 4 (Oar DO) = 20% EDO (D-6) 


In order to implement this derivation in a numerical 


analysis, equation (D-6) was rewritten as 


eS 2Q, + r,DQ, DS) 


Performing the change of independent variable (Appendix C) 


to accommodate a nonuniform mesh 


* —] * 
DQ; = (D xr.) “D Q(n,) (D-10) 


Substituting equations (D-8), (D-9) and (D-10) into 
SPauatcion (D-7) gives 


a 


— x 
= 2Q(n;) - r(n,)D_-x,) D Q(n,) (p=) 


a 


For the axis offset nonuniform mesh, r(n) is given by 
* . e 
equation (C-32) and (Dr) by equarronm(ea toe DUDScCLeueEend 


into equation (D-11) using equation (C-35) results in 
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2 


tanh[A(1-n,)] cosh [A(1=n,)] 7 
uw, = 20(ny) +1 - —poee yt i "in 
tanh[A(1-n,)] tanhicosh*[(1-n,)] . 
= I a Ss Q(n;) 


(een) 


For the wall offset mesh, equation (C-40) is substituted 
for equation (C~32) and all occurrences of the term l-n, 
are replaced by the term N° 


The value of u at the axis (uo) and at the wall 


(U4 7) 


were solved for by using the boundary conditions specified 


in Ref. 9, namely 


O(l) = 0 (Da 3) 
Donk) | ©=— A0 (D~14) 
Borg), y= Av (D= 1) 
376 K0) pee (D-16) 
From equations (D-13) and (D-14), using equation (D-7) 
memos oovious that 
u = 0 Cas) 


and from equations (D-15) and (D-7), it is Sima bar Ly 


found that 
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u =—e2zOr0y -; (D-18) 


where the finite difference approximation for Q(0) 1S given 


meeaquation (B-15). 
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